#----------------------------------------------------------------------------------------------------#
# Code for Main Text Analyses
# Engelhardt "Generational Persistence in the Nature of White Racial Attitudes"
#----------------------------------------------------------------------------------------------------#
library(stargazer)
library(lavaan)
library(xtable)
library(semTools)

# set working directory for folder with data and figures
# setwd("")

# Kam and Burge Analyses --------------------------------------------------
# Loading Data
kb_dat <- read.csv("./data/kam_burge.csv")


# * Traits of Blacks ------------------------------------------------------
RHS <- c("rr_sc*mil")

# Negative
m_negblk <- glm(as.formula(paste0("as.factor(neg_tblk)~", RHS)),
                data = kb_dat, family = binomial(link = "probit"))
summary(m_negblk)

# Positive
m_posblk <- glm(as.formula(paste0("as.factor(pos_tblk)~", RHS)),
                data = kb_dat, family = binomial(link = "probit"))
summary(m_posblk)


# * Individualism ---------------------------------------------------------
RHS <- c("rr_sc*mil")

# Affirmed
m_indvrule <- glm(as.formula(paste0("as.factor(indiv_rule)~", RHS)),
                  data = kb_dat, family = binomial(link = "probit"))
summary(m_indvrule)

# Flouted
m_indvflout <- glm(as.formula(paste0("as.factor(indiv_flouted)~", RHS)),
                   data = kb_dat, family = binomial(link = "probit"))
summary(m_indvflout)

# Broken
m_indvnotengh <- glm(as.formula(paste0("as.factor(indiv_notenough)~", RHS)),
                     data = kb_dat, family = binomial(link = "probit"))
summary(m_indvnotengh)

# * Discrimination ---------------------------------------------------------
RHS <- c("rr_sc*mil")

# Exits
m_discexist <- glm(as.formula(paste0("as.factor(discrim_prob_blk)~", RHS)),
                   data = kb_dat, family = binomial(link = "probit"))
summary(m_discexist)

# Denial
m_discnot <- glm(as.formula(paste0("as.factor(discrim_not_prob)~", RHS)),
                 data = kb_dat, family = binomial(link = "probit"))
summary(m_discnot)

# Reverse
m_discrev <- glm(as.formula(paste0("as.factor(discrim_reverse)~", RHS)),
                 data = kb_dat, family = binomial(link = "probit"))
summary(m_discrev)


# * Figure ---------------------------------------------------------------
# Sequence for X in predicted Values
p_rr <- seq(min(kb_dat$rr_sc, na.rm = T),
            max(kb_dat$rr_sc, na.rm = T),
            by = .1)
# Adjust
SHIFT <- .03


VARS <- c("neg_tblk", "pos_tblk", "indiv_rule", "indiv_flouted", "indiv_notenough", 
          "discrim_prob_blk", "discrim_not_prob", "discrim_reverse",
          "rr_sc", "mil")

# Predicted data
pred_dat_mill <- kb_dat[1:length(p_rr), VARS]
pred_dat_nmill <- kb_dat[1:length(p_rr), VARS]

# Millenials
for(i in 1:length(p_rr)){
  pred_dat_mill$neg_tblk[i] <- NA
  pred_dat_mill$pos_tblk[i] <- NA
  pred_dat_mill$indiv_rule[i] <- NA
  pred_dat_mill$indiv_flouted[i] <- NA
  pred_dat_mill$indiv_notenough[i] <- NA
  pred_dat_mill$discrim_prob_blk[i] <- NA
  pred_dat_mill$discrim_not_prob[i] <- NA
  pred_dat_mill$discrim_reverse[i] <- NA
  pred_dat_mill$rr_sc[i] <- p_rr[i]
  pred_dat_mill$mil[i] <- 1
}

# Non-Millenials
for(i in 1:length(p_rr)){
  pred_dat_nmill$neg_tblk[i] <- NA
  pred_dat_nmill$pos_tblk[i] <- NA
  pred_dat_nmill$indiv_rule[i] <- NA
  pred_dat_nmill$indiv_flouted[i] <- NA
  pred_dat_nmill$indiv_notenough[i] <- NA
  pred_dat_nmill$discrim_prob_blk[i] <- NA
  pred_dat_nmill$discrim_not_prob[i] <- NA
  pred_dat_nmill$discrim_reverse[i] <- NA
  pred_dat_nmill$rr_sc[i] <- p_rr[i]
  pred_dat_nmill$mil[i] <- 0
}

# Formating as Numeric
pred_dat_mill <- apply(pred_dat_mill, 2, as.numeric)
pred_dat_mill <- as.data.frame(pred_dat_mill)
pred_dat_nmill <- apply(pred_dat_nmill, 2, as.numeric)
pred_dat_nmill <- as.data.frame(pred_dat_nmill)


# ** Traits of Blacks ----------------------------------------
pred_negblk_mil <- predict(m_negblk, pred_dat_mill, se.fit=TRUE, type = "response")
pred_negblk_nmil <- predict(m_negblk, pred_dat_nmill, se.fit=TRUE, type = "response")
pred_posblk_mil <- predict(m_posblk, pred_dat_mill, se.fit=TRUE, type = "response")
pred_posblk_nmil <- predict(m_posblk, pred_dat_nmill, se.fit=TRUE, type = "response")


# pdf("./figures/pred_traits2_nocovs.pdf", width = 8, height = 3)
par(mfcol = c(1,3), mar = c(5,5,4,2))
plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Negative Trait Mention",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr, 
         y0 = pred_negblk_mil$fit - 1.96*pred_negblk_mil$se.fit, 
         y1 = pred_negblk_mil$fit + 1.96*pred_negblk_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_negblk_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT, 
         y0 = pred_negblk_nmil$fit - 1.96*pred_negblk_nmil$se.fit, 
         y1 = pred_negblk_nmil$fit + 1.96*pred_negblk_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_negblk_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1); axis(2, las = 1)

plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Positive Trait Mention",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr, 
         y0 = pred_posblk_mil$fit - 1.96*pred_posblk_mil$se.fit, 
         y1 = pred_posblk_mil$fit + 1.96*pred_posblk_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_posblk_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT, 
         y0 = pred_posblk_nmil$fit - 1.96*pred_posblk_nmil$se.fit, 
         y1 = pred_posblk_nmil$fit + 1.96*pred_posblk_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_posblk_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1, cex.axis = 1.25); axis(2, las = 1, cex.axis = 1.25)

plot(0:1, 0:1, type = "n",
     xlab = "",
     ylab = "",
     main = "",
     axes = F, cex.lab = 1.5)
# dev.off()


# ** Individualism -------------------------------------------
pred_indvrul_mil <- predict(m_indvrule, pred_dat_mill, se.fit=TRUE, type = "response")
pred_indvrul_nmil <- predict(m_indvrule, pred_dat_nmill, se.fit=TRUE, type = "response")
pred_indvflout_mil <- predict(m_indvflout, pred_dat_mill, se.fit=TRUE, type = "response")
pred_indvflout_nmil <- predict(m_indvflout, pred_dat_nmill, se.fit=TRUE, type = "response")
pred_indvnot_mil <- predict(m_indvnotengh, pred_dat_mill, se.fit=TRUE, type = "response")
pred_indvnot_nmil <- predict(m_indvnotengh, pred_dat_nmill, se.fit=TRUE, type = "response")

# pdf("./figures/pred_indiv2_nocovs.pdf", width = 8, height = 3)
par(mfcol = c(1, 3), mar = c(5,5,4,2))
plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Individualism Affirmed",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_indvrul_mil$fit - 1.96*pred_indvrul_mil$se.fit,
         y1 = pred_indvrul_mil$fit + 1.96*pred_indvrul_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_indvrul_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_indvrul_nmil$fit - 1.96*pred_indvrul_nmil$se.fit,
         y1 = pred_indvrul_nmil$fit + 1.96*pred_indvrul_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_indvrul_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1); axis(2, las = 1)

plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Individualism Flouted",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_indvflout_mil$fit - 1.96*pred_indvflout_mil$se.fit,
         y1 = pred_indvflout_mil$fit + 1.96*pred_indvflout_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_indvflout_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_indvflout_nmil$fit - 1.96*pred_indvflout_nmil$se.fit,
         y1 = pred_indvflout_nmil$fit + 1.96*pred_indvflout_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_indvflout_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1, cex.axis = 1.25); axis(2, las = 1, cex.axis = 1.25)

plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Individualism Broken",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_indvnot_mil$fit - 1.96*pred_indvnot_mil$se.fit,
         y1 = pred_indvnot_mil$fit + 1.96*pred_indvnot_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_indvnot_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_indvnot_nmil$fit - 1.96*pred_indvnot_nmil$se.fit,
         y1 = pred_indvnot_nmil$fit + 1.96*pred_indvnot_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_indvnot_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1, cex.axis = 1.25); axis(2, las = 1, cex.axis = 1.25)
# dev.off()


# ** Discrimination -------------------------------------------------------
pred_discexist_mil <- predict(m_discexist, pred_dat_mill, se.fit=TRUE, type = "response")
pred_discexist_nmil <- predict(m_discexist, pred_dat_nmill, se.fit=TRUE, type = "response")
pred_discnot_mil <- predict(m_discnot, pred_dat_mill, se.fit=TRUE, type = "response")
pred_discnot_nmil <- predict(m_discnot, pred_dat_nmill, se.fit=TRUE, type = "response")
pred_discrev_mil <- predict(m_discrev, pred_dat_mill, se.fit=TRUE, type = "response")
pred_discrev_nmil <- predict(m_discrev, pred_dat_nmill, se.fit=TRUE, type = "response")

# pdf("./figures/pred_disc2_nocovs.pdf", width = 8, height = 3)
par(mfcol = c(1, 3), mar = c(5,5,4,2))
plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Discrimination Exists",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_discexist_mil$fit - 1.96*pred_discexist_mil$se.fit,
         y1 = pred_discexist_mil$fit + 1.96*pred_discexist_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_discexist_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_discexist_nmil$fit - 1.96*pred_discexist_nmil$se.fit,
         y1 = pred_discexist_nmil$fit + 1.96*pred_discexist_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_discexist_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1); axis(2, las = 1)

plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Discrimination Denial",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_discnot_mil$fit - 1.96*pred_discnot_mil$se.fit,
         y1 = pred_discnot_mil$fit + 1.96*pred_discnot_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_discnot_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_discnot_nmil$fit - 1.96*pred_discnot_nmil$se.fit,
         y1 = pred_discnot_nmil$fit + 1.96*pred_discnot_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_discnot_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1, cex.axis = 1.25); axis(2, las = 1, cex.axis = 1.25)

plot(0:1, 0:1, type = "n",
     xlab = "Racial Resentment",
     ylab = "Predicted Probability",
     main = "Reverse Discrimination",
     axes = F, cex.lab = 1.5)
segments(x0 = p_rr,
         y0 = pred_discrev_mil$fit - 1.96*pred_discrev_mil$se.fit,
         y1 = pred_discrev_mil$fit + 1.96*pred_discrev_mil$se.fit,
         col = "black", lwd = 1.5)
points(p_rr, pred_discrev_mil$fit, col = "black", lwd = 1.5, pch = 15)
segments(x0 = p_rr+SHIFT,
         y0 = pred_discrev_nmil$fit - 1.96*pred_discrev_nmil$se.fit,
         y1 = pred_discrev_nmil$fit + 1.96*pred_discrev_nmil$se.fit,
         col = "grey", lwd = 1.5)
points(p_rr+SHIFT, pred_discrev_nmil$fit, col = "grey", lwd = 1.5, pch = 16)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 1)]), col = "black", side = 3)
rug(jitter(kb_dat$rr_sc[which(kb_dat$mil == 0)]), col = "grey")
legend("top", horiz = T, c("Millennials", "Older"), 
       col = c("black", "grey"), 
       lty = c(1,1), pch = c(15, 16), lwd = c(1.25,1.25), bty = "n")
axis(1, cex.axis = 1.25); axis(2, las = 1, cex.axis = 1.25)
# dev.off()

# * Table of Results (Appendix) ------------------------------------------------------
stargazer(m_negblk, m_posblk, m_indvrule, m_indvflout, m_indvnotengh, m_discexist, m_discnot, m_discrev,
          title = "Racial Resentment and Cognitive Responses to Items",
          model.numbers = F, 
          covariate.labels = c("Racial Resentment",
                               "Millennial",
                               "Racial Resentment*Millennial"),
          dep.var.caption = "",
          dep.var.labels = c("Negative Traits of Blacks", "Positive Trait of Blacks", 
                             "Individualism Affirmed", "Individualism Flouted", "Individualism Broken",
                             "Discrimination Exists", "Discrimination Denial", "Reverse Discrimination"),
          no.space = T, 
          notes.append = T,
          notes = c("Probit regression results. Standard errors in parentheses."), 
          notes.align = "l",
          out = c("./figures/KB_tab_nocovs.tex"),
          table.placement = "t",
          digits = 3, omit.stat = c("f", "adj.rsq"), align = T, df = F,
          label = "", header= F
)


# MG-CFA ------------------------------------------------------------------
# Loading Helper functions
source("./code/functions/SDI2.UDI2.RFunction.R")
source("./code/functions/invar_tab_function.R")

# Loading Data
anes16 <- read.csv("./data/anes16.csv")

# * Model -----------------------------------------------------------------
invar_dat <- subset(anes16, 
                    f2f == 1, 
                    select = c("specfavr", "pdisc", "thard", "dless", 
                               "millen"))

# Desired Fit Measures
FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

GROUP <- "millen"
mod_null <- '
# fixes intercepts across groups
specfavr ~~ c(p1, p1)*specfavr
pdisc ~~ c(p2, p2)*pdisc
thard ~~ c(p3, p3)*thard
dless ~~ c(p4, p4)*dless

# fixes variances across groups
specfavr ~ c(t1, t1)*1 
pdisc ~ c(t2, t2)*1 
thard ~ c(t3, t3)*1 
dless ~ c(t4, t4)*1 
'

fit_null <- cfa(mod_null, invar_dat, mimic = "Mplus", group = GROUP)
summary(fit_null)
null_fit <- fitMeasures(fit_null)[c("chisq","df", "cfi", "srmr", "rmsea")]


# Defining invariance model
mod <- '
rr =~ dless + thard + specfavr + pdisc 

# reverse wording
thard ~~ specfavr 
'

fit_config <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP)
summary(fit_config, standardized = T)
round(fitMeasures(fit_config)[FIT], 3)
config_fit <- fitMeasures(fit_config)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_metric <- cfa(mod, invar_dat, mimic = "Mplus",
                  group = GROUP, group.equal = c("loadings"))
summary(fit_metric, standardized = T)
round(fitMeasures(fit_metric)[FIT], 3)
metric_fit <- fitMeasures(fit_metric)[c("chisq","df", "cfi", "srmr", "rmsea")]

fit_scalar <- cfa(mod, invar_dat, mimic = "Mplus", 
                  group = GROUP,
                  group.equal = c("loadings", "intercepts"))
summary(fit_scalar, standardized = T)
round(fitMeasures(fit_scalar)[FIT], 3)
scalar_fit <- fitMeasures(fit_scalar)[c("chisq","df", "cfi", "srmr", "rmsea")]

comb_fit <- rbind(config_fit, metric_fit, scalar_fit)
comb_fit

## Permutation Method for significance
## fit indices of interest for multiparameter omnibus test
fit_measures <- c("chisq", "cfi", "srmr", "rmsea")

nperms <- 2000

## test configural invariance
set.seed(1693)
out_config <- permuteMeasEq(nPermute = nperms, 
                            con = fit_config, AFIs = fit_measures)
out_config

## test metric equivalence
set.seed(1693) # same permutations
out_metric <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_config, con = fit_metric, null = fit_null,
                            param = "loadings", AFIs = fit_measures)
summary(out_metric)

## test scalar equivalence
set.seed(1693) # same permutations
out_scalar <- permuteMeasEq(nPermute = nperms, 
                            uncon = fit_metric, con = fit_scalar, null = fit_null,
                            param = "intercepts", AFIs = fit_measures)
summary(out_scalar)

# * Substantive Differences ----------------------------------------------
MODEL <- fit_config
pars <- parameterEstimates(MODEL)
SDI2.UDI2(4, 
          pars$est[which(pars$op == "=~" & pars$group == 1)],
          pars$est[which(pars$op == "=~" & pars$group == 2)],
          pars$est[which(pars$op == "~1" & pars$lhs != "rr" & pars$group == 1)],
          pars$est[which(pars$op == "~1" & pars$lhs != "rr" & pars$group == 2)],
          sqrt(MODEL@Data@ov$var),
          pars$est[which(pars$op == "~1" & pars$lhs == "rr" & pars$group == 2)],
          sqrt(pars$est[which(pars$op == "~~" & pars$lhs == "rr" & pars$group == 2)]))

# * Tables -----------------------------------------------------------------
# Table 1
permute_fit <- array(NA, c(3, 12))
permute_fit <- as.data.frame(permute_fit)
rownames(permute_fit) <- c("Configural", "Metric", "Scalar")
colnames(permute_fit) <- c("$\\chi^2$", "CFI", "SRMR", "RMSEA", 
                           "$\\Delta\\chi^2$", "$\\chi^2$ p-value", 
                           "$\\Delta$CFI", "CFI p-value", 
                           "$\\Delta$SRMR", "SRMR p-value", 
                           "$\\Delta$RMSEA", "RMSEA p-value")
permute_fit["Configural", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["config_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["metric_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Metric", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_metric@AFI.obs
permute_fit["Metric", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_metric@AFI.pval

permute_fit["Scalar", 
            c("$\\chi^2$", "CFI", "SRMR", "RMSEA")] <- comb_fit["scalar_fit", c("chisq", "cfi", "srmr", "rmsea")]
permute_fit["Scalar", 
            c("$\\Delta\\chi^2$", "$\\Delta$CFI", "$\\Delta$SRMR", "$\\Delta$RMSEA")] <- out_scalar@AFI.obs
permute_fit["Scalar", 
            c("$\\chi^2$ p-value", "CFI p-value", "SRMR p-value", "RMSEA p-value")] <- out_scalar@AFI.pval

tab <- xtable(permute_fit, digits = 3, 
              caption = "Measurement Invariance of Racial Resentment, Millennials vs. Older Whites",
              align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c", "c"),
              label = "anes_invar",
              # specifying significant digits
              display = c("s", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg", "fg"))
NOTE <- list()
NOTE$pos <- list()
NOTE$pos[[1]] <- c(nrow(tab))
NOTE$command <- c("\\bottomrule \\multicolumn{10}{l}{\\multirow{2}{\\textwidth}{Note: The configural model freely estimates item loadings using the \\textit{deserve less} item to define the dimensions. The metric model constrains item loadings to equality. The scalar model constrains both item loadings and intercepts to equality. Data from face-to-face interviews in the 2016 ANES.}}")
print(tab, caption.placement = "top",
      sanitize.text.function = function(str) gsub("$\backslash$", "\\", str, fixed = TRUE),
      add.to.row = NOTE,
      hline.after = c(-1, 0),
      file = "./figures/anes_invar.tex")

# Table 2
invar_tab(config = fit_config, 
          metric = fit_metric, 
          scalar = fit_scalar, 
          title = "Generation Invariance, Face-to-Face Interviews ANES 2016",
          file_name = "./figures/anes_f2f_models_tab.tex", 
          varnames = c("Deserve Less", "Try Hard", "Special Favors", "Past Discrimination"))

